
options(rasterTmpDir='C:/temp/')

##
## IMPORTANT PATH to Rtools to run Cpp
##
Sys.setenv(PATH = paste("C://WBG//Utility//Rtools//bin", Sys.getenv("PATH"), sep=";"))
Sys.setenv(BINPREF = "C://WBG//Utility//Rtools//mingw_$(WIN)//bin//")

require(Rcpp)

output_version <- '2020-10-20'

ntl.zones <- read_sf(dsn=paste0(wkdir,'/local/gis_data/003_boundaries/fishnet'),layer='fishnet_lake_chad_ntl_adm0_x_d01')

ntl.zones.prj <- ntl.zones %>%
  # 102022
  st_transform(crs="+proj=aea +lat_1=20 +lat_2=-23 +lat_0=0 +lon_0=25 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m no_defs") %>% 
  st_centroid() %>% 
  st_geometry()

M.ntl.zones <- ntl.zones.prj %>% 
  st_transform(4326) %>%
  st_coordinates() 
head(M.ntl.zones)

# convert polygons to points 
# need to check border cases -> INSIDE or NOT? #
ntl.zones.pts = st_as_sf(as.data.frame(M.ntl.zones), coords = c("X", "Y"), crs = 4326)
ntl.zones.tb = st_drop_geometry(subset(ntl.zones,select=c(OBJECTID, GID_0)))
ntl.zones.pts = dplyr::bind_cols(ntl.zones.pts,ntl.zones.tb)

##
## BORDER POST : JEDWAB AND STOREYGARD 2020
## 2020-04-20
##
bpost <- st_read(paste0(wkdir,'/local/gis_data/018_transportation/border_crossings_chad_basin_2008_v1.shp'))

## 
bpostgid0 <- st_join(bpost,subset(ntl.zones,select=c(OBJECTID,GID_0)))

## correct assignment
ntl.zones.pts$dist2bcrossing<- NA

# CMR 
ntl.zones.pts$dist2bcrossing[ntl.zones.pts$GID_0=="CMR"] <- unlist(nngeo::st_nn(ntl.zones.pts[ntl.zones.pts$GID_0=="CMR",], subset(bpostgid0,GID_0=="CMR"), k = 1, returnDist = T)$dist) / 1000

# NER
ntl.zones.pts$dist2bcrossing[ntl.zones.pts$GID_0=="NER"] <- unlist(nngeo::st_nn(ntl.zones.pts[ntl.zones.pts$GID_0=="NER",], subset(bpostgid0,GID_0=="NER"), k = 1, returnDist = T)$dist) / 1000

# NGA
ntl.zones.pts$dist2bcrossing[ntl.zones.pts$GID_0=="NGA"] <- unlist(nngeo::st_nn(ntl.zones.pts[ntl.zones.pts$GID_0=="NGA",], subset(bpostgid0,GID_0=="NGA"), k = 1, returnDist = T)$dist) / 1000

# TCD
ntl.zones.pts$dist2bcrossing[ntl.zones.pts$GID_0=="TCD"] <- unlist(nngeo::st_nn(ntl.zones.pts[ntl.zones.pts$GID_0=="TCD",], subset(bpostgid0,GID_0=="TCD"), k = 1, returnDist = T)$dist) / 1000

save.dta13(st_drop_geometry(subset(ntl.zones.pts,select=c(OBJECTID,GID_0,dist2bcrossing))), 
           paste0(wkdir,"/proc_data/DIST2BCROSS",output_version,".dta"))